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Abstract. We report a molecular dynamics simulation of melting of tungsten (W) nanoparticles. The 
modified embedded atom method (MEAM) interatomic potentials are used to describe the interaction 
between tungsten atoms. The melting temperature of unsupported tungsten nanoparticles of different 
sizes are found to decrease as the size of the particles decreases. The melting temperature obtained in 
the present study is approximately a decreasing function of inverse radius, in a good agreement with the 
predictions of thermodynamic models. We also observed that the melting of a W nanoparticle is preceded 
by the premelting of its outer skin at a temperature lower than its melting temperature. 

PACS. 61.46.Df Nanoparticles - 61.46.-w Nanoscale materials - 64.70.Dv Solid-liquid transitions - 
65.80.-f n Thermal properties of small particles, nanocrystals, and nanotubes 



' 1 Introduction 

. Tungsten, along with its alloys and compounds, occupies a 

■ unique position in materials science. The material proper- 
] ties that make tungsten attractive to the metals industry 

■ are high density, hardness, melting temperature, elastic 
modulus, and conductivity in conjunction with the low 
thermal expansion. The combination of these unique prop- 

■ erties explains the diverse applications of tungsten rang- 
ing from home lighting to thermonuclear fusion first-wall 
protection [1,2]. With nanoscale tungsten powders avail- 
able at a reasonable cost, its usage will increase greatly 

, and a new approach is required to balance the size de- 
' pendent advantages against the temperature dependent 
limitations. Therefore, it is of great importance to un- 
1 derstand the thermal stability of tungsten nanoparticles 
for their applications at higher temperatures. It has been 
seen previously that nanoparticles exhibit significant de- 
crease in melting temperatures compared to infinite bulk 
solids [3]. This is related to the fact that the liquid- vapor 
interface energy is generally lower than the average solid- 
vapor interface energy [4]. Based on thermodynamics, a 
phenomenological relation between the melting tempera- 
ture and particle size has been obtained: the melting tem- 
perature of a nanoparticle decreases as inversely propor- 
tional to the particle diameter [3-5] . It is also known that 
premelting, the phenomenon where the surface atoms of 
nanoparticles lose their solid ordering and hence melt prior 
to complete melting of the whole particle [5-10], plays an 
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important role in understanding the melting of nanoparti- 
cles. On the other hand, the Hall-Pctch effect — the hard- 
ness increases in proportional to the inverse square-root 
of grain size [11,12] — suggests that significant opportuni- 
ties exist if nanoscale powders could be consolidated to 
full densities with minimized coarsening. Hence, knowl- 
edge of accurate melting temperature for different parti- 
cle size plays an important role for the advancement of 
present engineering and technological growth. 

Molecular dynamics simulations offer an effective tool 
to study the melting and coalescence of nanoparticles [13, 
14]. These atomistic simulations require accurate atomic 
interaction potentials to compute the total energy of the 
system. First-principles calculations can provide the most 
reliable interatomic potentials. However, realistic simula- 
tions of the melting of nanoparticles often require a num- 
ber of atoms that renders these methods impractical: they 
either require too much computer memory or take too long 
to be completed in a reasonable amount of time. One al- 
ternative is to use empirical or semi-empirical interaction 
potentials that can be evaluated efficiently. In this study, 
we use the modified embedded atom method (MEAM) 
originally proposed by Baskes et al. [15, 16]. MEAM was 
the first semi-empirical atomic potential using a single 
formalism for fee, bcc, hep, diamond-structured materi- 
als and even gaseous elements, in good agreement with 
experiments or first-principles calculations [16, 17]. The 
MEAM is an extension of the embedded-atom method 
(EAM) [18-20] to include angular forces. Cherne et al. 
made a careful comparison of MEAM and EAM calcula- 
tions in a liquid nickel system [21]. 
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Atomistic simulations of a wide range of elements and 
alloys have been performed using the MEAM potentials. 
A realistic shear behavior for silicon was first obtained us- 
ing the MEAM by Baskcs ct al. [15]. The MEAM was also 
applied to various single elements [16] and to silicon-nickel 
alloys and interfaces [22]. Gall et al [23] used the MEAM to 
model the tensile debonding of an aluminum-silicon inter- 
face. Lee and Baskes [24] extended the MEAM to include 
the second nearest-neighbor interactions. A new analytic 
modified embedded-atom method (AMEAM) many-body 
potential was also proposed and applied to several hep 
metals, including Mg [25,26]. For the Mg-Al alloy sys- 
tem, a set of EAM potentials has been developed using 
the "force matching" method by Liu ct al [27]. Recently, a 
new set of MEAM potentials for Mg-Al alloy system was 
developed by Jelinek et al [28]. These new potentials show 
a significant improvcmcnit over the previously published 
potentials, especially for the surface formation, stacking 
faults, and point defect formation energies. 

The paper is organized in the following manner. In 
Sec. 2, we give a brief review of the MEAM. In Sec. 3, the 
procedure for melting simulation is presented. MD simula- 
tion results are presented and discussed in Sec. 4. Finally, 
in Sec. 5, we summarize our findings. 



2 MEAM theory 

The total energy E oia system of atoms in the MEAM [29] 
is approximated as the sum of the atomic energies 



E : 



(1) 



The energy of atom i consists of the embedding energy 
and the pair potential terms: 



i^i =i^^(p.) + ^E^»J"('^»^")• 



(2) 



Fi is the embedding function of atom i, pi is the back- 
ground electron density at the site of atom «, and (t>ij{rij) 
is the pair potential between atoms i and j separated by a 
distance rjj. The embedding energy Fi(pi) represents the 
energy cost to insert atom i at a site where the background 
electron density is pi. The embedding energy is given in 
the form 

F,Xp,) = A,E^f5,\n{pi). (3) 

where the sublimation energy and parameter Aj de- 
pend on the element type of atom i. The background elec- 
tron density pi is given by 
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and 
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The zcroth and higher order densities, p^f^ , pt\ pf^ ■< 
Pi are given in Eq. (9). The composition-dependent elec- 
tron density scaling p^ is given by 



9° = p,oZ,oG{rf), 



(7) 



where pio is an element-dependent density scaling, 

is the first nearest-neighbor coordination of the reference 

system, and /Y**^ is given by 



^k) (k) 
°i 1 



(8) 



fc=i 



(k) 

where s\ is the shape factor that depends on the refer- 
ence structure for atom i. Shape factors for various struc- 
tures are specified in the work of Baskes [16]. The partial 
electron densities are given by 
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where 



I] a 



is the a. component of the displacement vec- 



tor from atom i to atom j. Sij is the screening function 
between atoms i and j and is defined in Eqs. (16). The 
atomic electron densities are computed as 
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where r^* is the nearest-neighbor distance in the single- 

(k) 

element reference structiirc and is element-dependent 
parameter. Finally, the average weighting factors are given 

by 



Pt 



(11) 



where tj*^^ is an element-dependent parameter. 
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The pair potential is given by 



(12) 
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(15) 

where tty is an element-dependent parameter. The subli- 
mation energy . the equilibrium nearest- neighbor dis- 
tance r? , and the number of nearest-neighbors Zij are 
obtained from the reference structure. 

The screening function Sij is designed so that Sij = 
1 if atoms i and j are unscreened and within the cutoff 
radius rc, and Sij = if they are completely screened or 
outside the cutoff radius. It varies smoothly between and 
1 for partial screening. The total screening function is the 
product of a radial cutoff function and three-body terms 
involving all other atoms in the system: 
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Note that Cmin and Cmax can be defined separately for 
each i-j-k triplet, based on their element types. The pa- 
rameter Ar controls the distance over which the radial 
cutoff function changes from 1 to near r = rc- 



3 Molecular dynamics simulation 
3.1 Atomic potential 

We use the MEAM potential parameters for tungsten (W) 
proposed by Baskes [16]. The potential parameters that 
are used for our simulation of W nanoparticles are listed 
in Table 1. These parameters are obtained by fitting the 
room temperature elastic properties using bcc as the ref- 
erence structure. Cmax and Cmin are chosen to consider 
only the first nearest-neighbor interactions [30]. 



We validate the potential by computing different phys- 
ical properties of tungsten systems and comparing them 
with DFT calculations. The results are compared with 
those of DFT calculations as shown in Table 2. Energy 
calculations and geometry optimizations of various struc- 
tures were performed using Blochl's all-electron projector 
augmented wave (PAW) method [31] as implemented by 
Kresse and Joubert [32]. For the treatment of electron 
exchange and correlation, we use the generalized gradi- 
ent approximation (GGA) using Perdew-Burke-Ernzerhof 
scheme [33]. 



3.2 Simulation Procedure 

We performed a detailed MD simulation of the melting of 

unsupported spherical bcc W nanoparticles, 2 12 nm in 
diameter (259-56905 atoms). The surface boundary con- 
dition was free and no external pressure was applied. Each 
nanoparticle was constructed by cutting out atoms within 
a specified radius from the tungsten bulk in bcc structure. 
The equations of motion were integrated using time steps 
At = 4 X 10~^^ s. We begin each MD run by random- 
izing the atomic velocities of the nanoparticle according 
to the Maxwell-Boltzmann distribution. We increase the 
temperature of the heat bath in steps of AT — 100 K from 
the initial temperature Ti = 500 K to the final tempera- 
ture up to Tf = 4000 K. We let the nanoparticles come to 
equilibration for 50 000 time steps at each temperature. 
Statistical (time-averaged) data for the energetics are col- 
lected after the system has adjusted to the new tempera- 
ture, which is typically after 25 000 time steps following 
a temperature increase. For the particles of diameters less 
than 8 nm, 20 000 time steps were used to adjust the par- 
ticles to each new temperature. The isothermal condition 
was maintained by using Nose-Hoover thermostat [34,35]. 



4 Results and Discussion 

The most straightforward method to identify the melt- 
ing of atomistic structures in MD simulations is to moni- 
tor the variation of the internal energy with temperature. 
Fig. 1 shows the internal energies of the W nanoparti- 
cles with different diameters as a function of temperature. 
It is clearly seen from the Fig. 1 that each internal en- 
ergy curve goes from one linear region to another. The 
overall melting is clearly identified by the abrupt "jump" 
in the internal energy curve, the height of the jump is a 
measure of AHm, the amount of heat required for melt- 
ing, and it decreases as the size of nanoparticle decreases. 
The melting temperatures calculated based on the present 
MD simulation of W nanoparticles are listed in Table 3. 
We note that the melting temperature of bulk W from 
our MD simulation, 3900 K, is slightly higher than the 
experimentally measured value of 3695 K [36]. The dis- 
crepancy in this result is mainly due to the super-heating 
of the simulated lattice, as it has been observed that the 
confined lattice without free surface can be significantly 
superheated [37,38]. Although, it is not the main focus 
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of this study, one can follow the procedure prescribed by 
Morris et al [39] to establish co-existence of solid and liq- 
uid phases to determine the melting temperature of the 
bulk W system without supcr-hcating. More importantly: 
we also note that the melting temperature decreases dras- 
tically as the size of the particle decreases. This result 
suggests that the thermal stability of small nanoparticles 
must be carefully investigated before they can be used in 
applications such as nano-devices. 

The melting behavior of 2 nm particle seems to be 
different from those of bigger particles: at the onset of 
the melting, the internal energy curve dips down before 
climbing up again. A similar behavior has been observed 
in the melting of small Au nanoparticles of diameters up to 
2.8 nm [13]. The onset of melting provides surface atoms 
an opportunity to rearrange themselves to optimize the 
local morphology and lower their portion of the internal 
energy. For extremely small particles, where the surface 
area to volume ratio is large, this will cause the total in- 
ternal energy of the particle to decrease briefly as shown 
in Fig. 1. However, a further detailed study focusing on 
small nanoparticles will be required to elucidate this pe- 
culiar behavior. 

The variation of the melting temperature with the size 
of the W nanoparticles is plotted in Fig. 2. The melt- 
ing point depression of W nanoparticles exhibit the same 
qualitative behavior found in the MD simulation of Au 
nanoparticles [9,13]. A similar size dependence of melt- 
ing point depression has been observed experimentally 
over a broad range of particle sizes for particles in cluster 
beams [40-43] as well as particles on substrates [6,7,44,45]. 

For spherical particles of diameter R, a melting tem- 
perature Tm{R) can be obtained phenomenologically [4,5, 
46] by equating the Gibbs free energies of solid and liquid 
spherical clusters, assuming constant pressure conditions: 

r„(i?) = r4(i-|-), (17) 

where is the melting temperature of the bulk tung- 
sten and i?i is a parameter related to physical quanti- 
ties such as the solid and liquid densities, the bulk la- 
tent heat of melting, and solid-vapor and liquid- vapor in- 
terface energies. In obtaining this model, the surface en- 
ergy anisotropy of the solid is not taken into account, and 
the possibility of inhomogeneous phases (such as a liquid 
layer due to premelting) is also neglected. The solid line in 
Fig. 2 corresponds to the simple thermodynamical model 
of Eq. (17), with constant parameters T* = 2900 K and i?i 
= 1.5 nm. The curve shows clearly that the melting point 
of W nanoparticles decrease according to 1/i? dependence 
as predicted in Eq. (17). However, the fitted value of T* 
is significantly lower than the melting temperature of the 
bulk tungsten. This result indicates that the characteris- 
tics of the curve is likely to change for nanoparticles with 
larger diameters, and further study with larger nanopar- 
ticles will be needed to test the applicability of this model 
to W nanoparticles. 

Hanszen [47] proposed another model of melting in 
terms of classical thermodynamics assuming that a liq- 



uid over-layer forms at the surface of the solid cluster 
and grows towards the solid core, below the melting point 
[48,49]. When the liquid layer thickness exceeds a critical 
thickness, the whole cluster melts homogeneously. In this 
model, the melting point Tm{R) of W nanoparticles with 
diameter R can be expressed as [6, 7] 

T„(i?)=T* (l--^ + :|), (18) 

When the data of Table 3 were fitted to Eq. (18), we ob- 
tained negligibly small values for to and R2 , thus returning 
to the model of Eq. (17). 

Fig. 3 shows the cross sections of a W nanoparticle 
with a diameter 10 nm through the center of the parti- 
cle. The displacement vectors at different times during the 
MD simulation run at the temperature 2000 K are plotted. 
Fig. 3 shows that at a temperature below the melting point 
the atoms in the entire nanoparticle vibrate in their places 
while retaining their bee crystal structure. As the temper- 
ature increase, several layers of atoms start to lose their 
periodicity and form a liquid shell as shown in Fig. 3(b). 
Once the thickness of the liquid layer reaches a critical 
thickness, the whole nanoparticle melts. Our MD simula- 
tion confirms the experimental observation that nanoscale 
materials simultaneously display solid-like and liquid-like 
characteristics, and under thermodynamic equilibrium, a 
fraction of the atoms in the outer shell of the particle ex- 
hibit liquid- like behavior and the remaining fraction in the 
inner core act as solid [7]. Hence, melting point depression 
and the presence of disorder in nanoscale W powders will 
play an important role in various industries, including mi- 
croelectronic industries such as printed circuit board drill 
bits [50]. 

Fig. 4 shows a few snapshots of MD simulation for a 
small W nanoparticle with the diameter of 2 nm. We note 
that our simulations does not show pronounced faceted or 
step-like structures. We found a similar result when the 
nanoparticles are heated to the melting temperature and 
cooled down slowly. Our results are in good agreement 
with an earlier experiment that found no evidence for a 
faceted or step-like microstructure in a single tungsten 
crystal [51]. 



5 Conclusions 



The thermal stability of unsupported W nanoparticles has 
been investigated using a MD simulation. The MEAM po- 
tential was used to described the interatomic interactions. 
W nanoparticles melt at a temperature that is lower than 
the bulk melting temperature. The result of our present 
calculation shows that the melting temperature to be ap- 
proximately a decreasing function of inverse radius. We 
found that W nanoparticle melting is preceded by surface 
melting effects of its outer skin, similar to the melting of 
spherical clusters of many other elements. 
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Table 2. Calculated physical properties of W using the present 
MEAM parameters in comparison with DFT calculations. Bo 
is the bulk modulus (GPa); Cn, C12, C44 are the elastic con- 
stants (GPa); -E(ioo), -E{iio)) ^'(iii) are surface energies of cor- 
responding surfaces (mj/m^); AE's are the structural energy 
differences (eV/atom). 



Parameter DFT MEAM 



Bo 


330 


270 


(Cii - Ci2)/2 


190 


160 


C44 


280 


160 


-£■{100) 


7810 


5980 


-E{iio) 


6390 


5660 


£^{111) 


7190 


5030 


■^-^bcc — 'fee 


0.494 


0.325 


'^-S'bcc — *hcp 


0.397 


2.168 



Table 3. Melting temperatures of W nanoparticles with dif- 
ferent diameters 



Diameter (nm) 


No. of atoms 


(K) 


2.0 


259 


1000 


4.0 


2085 


1900 


6.0 


7119 


2200 


8.0 


16865 


2300 


10.0 


33079 


2400 


12.0 


56905 


2500 


Bulk 


00 


3900 




1000 2000 3000 4000 5000 



Temperature (K) 



Fig. 1. Internal energies of the W nanoparticles with different 
diameters as a function of temperature. The same data for W 
bulk are also shown. 
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Table 1. The MEAM potential parameters for W from Ref. 30. E is the sublimation energy, r*^ is the equilibrium nearest- 
neighbor distance, A is the scaling factor for the embedding energy, a is the exponential decay factor for the universal energy 
function, /j'""^) are the exponential decay factors for the atomic densities, t^^~^'^ are the weighting factors for the atomic 
densities, Cmax and Cmin are the screening parameters. 



E° [eV] 


r« [A] A 


a 






^(3) ^(0) ^(1) ^(2, 


+(3) 1^ n . 

f' ^max ^min 


8.66 


2.74 0.98 


5.63 3.98 


1.00 


1.00 


1.00 1.00 3.16 8.25 


-2.70 2.8 2.0 




Fig. 3. (color online). The cross section of W nanoparticle with 
a diameter 10 nm through the center of the particle showing the 
displacement vectors in the interval of 6 ps at the temperature 
of (a) 300 K and (b) 2000 K. The color and the size of the 
spheres represent the magnitude of the displacement vectors. 




Fig. 4. Snapshots of MD simulation at (a) 300 K, (b) 900 K, 
and (c) 1500 K (above the melting temperature) for 2 nm par- 
ticle. 



